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Q^; Abstract 



We present two open-source (BSD) implementations of ellipsoidal harmonic expansions for solv- 
ing problems of potential theory using separation of variables. Ellipsoidal harmonics are used 
surprisingly infrequently, considering their substantial value for problems ranging in scale from 



< 

u 

' molecules to the entire solar system. In this article, we suggest two possible reasons for the paucity 

relative to spherical harmonics. The first is essentially historical — ellipsoidal harmonics developed 

(N . 

^ , during the late 19th century and early 20th, when it was found that only the lowest-order harmon- 

ic 

' ics are expressible in closed form. Each higher-order term requires the solution of an eigenvalue 

(N ■ 

■ problem, and tedious manual computation seems to have discouraged applications and theoretical 

. studies. The second explanation is practical: even with modern computers and accurate eigenvalue 

(N 

algorithms, expansions in ellipsoidal harmonics are significantly more challenging to compute than 
those in Cartesian or spherical coordinates. The present implementations reduce the "barrier to 
5h ■ entry" by providing an easy and free way for the community to begin using ellipsoidal harmonics in 

actual research. We demonstrate our implementation using the specific and physiologically crucial 
problem of how charged proteins interact with their environment, and ask: what other analytical 
tools await re-discovery in an era of inexpensive computation? 
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I. INTRODUCTION 



It is intuitively obvious that a cow is much better approximated by an eUipsoid than by a 
sphere. Less obvious, or at least less well known, is the method of ellipsoidal harmonics for 
solving the ellipsoidal cow exactlyii^^. In an effort to increase the popularity and impact of 
these fascinating functions, we present in this paper two open-source (BSD) implementations 
for calculating the ellipsoidal harmonics and solving problems of potential theory (they may 
be downloaded freely online^). To demonstrate the methods' correctness, we model the 
physiologically crucial electrostatic interactions between a protein and surrounding water 
using a popular continuum theory based on the Poisson equation^'^. Our results indicate, 
that improved numerical methods are needed, and we hope that releasing these codes will 
not only encourage application scientists to try ellipsoidal harmonics, but also motivate 
numerical analysts to help develop more accurate computational approaches. 

Our implementations, written in Python and MATLAB^, employ a variety of insights 
developed over more than a century of research^'^'^; many important results were published 
during the mid- 19th century and early 20th by greats like Heine^, and Charles Darwin's son 
Sir George Howard Darwin^°. From the authors' perspective, these surprisingly direct ties 
to mathematical history offer an unusual and inspiring aspect to their study. Ellipsoidal 
coordinates are the most general coordinate system for which the Laplace equation is sepa- 
rable, as the other coordinate systems for which Laplace is separable can all be considered 
degenerate forms of ellipsoidal coordinates^-. In fact, the functions for each coordinate di- 
rection satisfy the same scalar equation, the Lame equation. It seems somewhat ironic that 
the most general separable coordinate system should lead to this simple Cartesian-like result 
rather than to a more complicated form as found e.g. in spherical harmonics; the rabbit 
hole goes much deeper, but for now it suffices to observe that ellipsoidal harmonics are not 
a simple generalization of spherical harmonics. 

Unsurprisingly, the more general shape allows ellipsoidal-harmonic expansions to be 
more accurate than ones based on spherical harmonics^'^^'i^. Successful applications cover 
most of potential theory: gravity-ii^, electrostatics^*^"— , electromagnetics^"—, hydrody- 
namics^^"— , and elasticity^. Specific uses in molecular and biological sciences include solv- 
ing the Schrodinger equation^^, modeling van der Waals (close packing) interactions between 
molecules^, design of magnetic resonance imaging (MRI) devices2§. and analysis of clinical 
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electroencephalography (EEG) and magnetoencephalography (MEG) data^ii^^. Unfortu- 
nately, despite the range of applications and their advantages over spheres, a comparatively 
small number of publications actually use, or encourage the use of, ellipsoidal harmonics in 
practice. Notable and welcome exceptions may be found in Sten's excellent tutorial review^^ 
and presumably in the upcoming book of Dassios^, who has been one of the subject's leading 
developers and championsi^i^ii^-J^. A few older books also address the theory^^*^, but for 
decades the authoritative reference has been E. W. Hobson's 1931 text^, published only two 
years before he passed away (for a charming obituary, see^). 

One reason for the relative dearth of papers using ellipsoidal harmonics is that they are 
hard to compute: unfortunately for investigators who might be interested in actually using 
ellipsoidal harmonics, widely available references do not address the numerous critical chal- 
lenges associated with actually computing ellipsoidal harmonics to arbitrary order. Naive 
numerical algorithms for their determination are unstable^^, and careful reformulations are 
requirecl^>i^i^~— . This difficulty contrasts sharply to calculations in spherical harmonics or 
Cartesian space, where the functions for the expansion are given in closed form and can be 
computed easily. A second possible explanation, more speculative, is that the foundations 
of ellipsoidal harmonic analysis were developed many decades before the advent of digital 
computersii^i^; actually using ellipsoidal harmonics to arbitrary order necessitates exten- 
sive computation, whether manual or by computer. Computation is necessary because the 
harmonics are defined as polynomials whose coefficients are solutions to eigenvalue problems. 
As such only the lowest-order modes can be expressed in closed form^, and all other modes 
must be determined numerically (by the Abel-Ruffini theorem). Furthermore, the harmonics 
depend on the semi-axes of the ellipsoid, necessitating a new set of tedious manual compu- 
tations for every problem. For that matter, we note that even the modes available in closed 
form require cumbersome algebraic manipulations^*^. It seems that the large quantities of 
manual arithmetic led practitioners in many fields to adopt less accurate, but much more 
easily computed, spherical harmonic methods. 

Thus, it seems possible that open-source implementations of ellipsoidal harmonics, based 
on the latest numerical algorithms for their computation, could encourage their re-discovery 
and wide application in computational science. The present work offers a simple demonstra- 
tion that deriving an analytical solution using ellipsoidal harmonics is not any harder than 
deriving one using spherical harmonics. Instead, the challenge is in computing the functions 
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themselves, a problem met by the released implementations; though, again, we acknowledge 
that a number of improvements are needed. We begin the technical portion of the paper in 
Section [III which presents the application setting (continuum electrostatics for biomolecular 
solvation) and the mathematical components needed to compute ellipsoidal harmonics. In 
Section UTTl we describe our implementations at a high-level but highlight critical details and 
subtle issues. Section HVl presents computational results, and Section |V] concludes the paper 
with a discussion of remaining challenges and possible future directions. 



II. THEORY 



Space constraints preclude full accounting of these topics; we present only the most 
relevant formulae and refer interested readers to the extensive and excellent treatments 
which have made our work possible^'^'^'^i. It is unfortunate that the clearest publications do 
not use a single, consistent notation; some types of analyses and identities are more readily 
seen in one notation or another. We have attempted to use the most popular conventions 
where possible, and apologize for introducing yet another mongrel notation. 



A. Biomolecule Electrostatics 



We assume that the biomolecule is a triaxial ellipsoid whose surface satisfies 

xp' 
0^ 

with 0<c<6<a<oo called the ellipsoid's semiaxes. In this paper, we use the term 
ellipsoid to mean a tri-axial ellipsoid (three unequal semi-axes), which is the most general 
ellipsoidal geometry. Ellipsoids with two equal semi-axes are called spheroids, these special 
cases are still more more straightforward to study^"— and apply^"— . 

The molecular interior is assumed to be a homogeneous local dielectric of relative per- 
mittivity ei obeying linear response, and the molecular charge distribution is (without loss 
of generality) assumed to be a set of Q discrete point charges, the ith of which being situ- 
ated at Yi and having value qi. Thus, inside the molecule (called region 1), the electrostatic 
potential obeys the Poisson equation 

1 ^ 

V2$i(r) = Vg,5(r-r,) (2) 

^1 k=i 



where S{r) is the Dirac delta function. The solvent outside the molecule (region 2) is modeled 
as a homogeneous dielectric with dielectric constant €2 and no fixed charges (i.e. we are 
modeling a salt-free solution), so the potential in this region obeys the Laplace equation 

V2$2(r) = 0. (3) 

At the dielectric boundary, i.e., the elhpsoid surface defined by Eq. [T], the potential is 
continuous, as is the normal component of the electric displacement field D(r) = e(r)E(r), 
so for a point on the surface, 

<!>,{rs) = Mrs)) (4) 

and we assume that the normal direction n(r5) points outward into the solvent. The poten- 
tial is assumed to decay sufficiently rapidly as ||r|| — )• 00. We note that more sophisticated 
electrostatic models have been presented, e.g.-; here we use a simple and popular one^ to 
demonstrate our implementation. 



B. Ellipsoidal Coordinates and Separation of Variables 

In the ellipsoidal coordinate system, a point denoted in Cartesian space r = {x, y, z) is 
written as (A,/i, z/); each ellipsoidal coordinate is a root of the cubic 

with 

h^ = a^- (7) 
k^ = a^- (8) 

and we take the positive square roots so that < h < k. The squares of the ellipsoidal 
coordinates are in the ranges 

A' e [e, +00 ) (9) 

/i' G [h^, k^] (10) 

G [0, h'] . (11) 



z/2 
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Points on the surface of the elhpsoid with semi-axes a, b, and c satisfy X = a. Many 
texts enforce that A is positive, i.e. A G [k,+oo ), by analogy with the radial coordinate 
in spherical systems. We have found non-negativity assumption problematic for inverse 
coordinate transforms and do not use it (see discussion in Sections IIIII and |V|) . 

Our expressions for coordinate transformations come from Romain and Jean-Pierre^, who 
note that for a given point r = (x, y, z), the magnitudes of the ellipsoidal coordinates (A, /i, z/) 
can be computed via 

A2 = 2v^cos('?V^ (12) 



At =2VQcos - + — 



9 r— ,'9 2tt\ Wi 
iy2 = 2VQcos -+ ^ 



3J 3 

6 Att\ 

3 ^ Ty 3 

e 27t\ Wi 

3 ^ Ty T 



where 



Q = "^^^^ (13) 
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9wiW2 - 27ws - 2wl 
R = ^4 (14) 



R 

cose = ^= (15) 



= -(x^ + ?/^ + ;2^ + /i2 + A;2) (16) 
yj, = x\h^ + k^) + y^e + z^h^ + h^e (17) 
= -x'h''k\ (18) 



The Cartesian coordinates can be computed from the ellipsoidal ones via 



AVV2 



X 



2 



y 



h^k^ 

~ h^{k^-h^) 
2 (A2-P)(P-^2)(p_^2 



(19) 



^ P(A;2 - h^) 

Romain and Jean-Pierre note that the simple Cartesian-to-ellipsoidal transformation of 
Eq. [12] suffers accuracy problems that can be important in special cases, and present a 
more sophisticated approach for improved accuracy^. Our present implementation includes 
only Eq. ^2\ and the improved algorithm is under development. All coordinate transfor- 
mations using this simpler method are verified by computing the inverse transform and 



checking with a tolerance of 10 ^. Also, we note that the above expressions correct a small 
typographical error in Eq. 7 of their work. 



C. The Lame Equation and its Solutions 

In ellipsoidal harmonics, the Laplace equation separates such that the solutions for each 
coordinate satisfy the same differential equation (Eq. [201) . which is called Lame's equation: 

(,2 _ i^2^^^2 _ k^)^^s) + s{2s^ -h'- e)^{s) + {p- qs')E{s) = (20) 
as^ as 

where p and q are unknown constants. It turns out that q = n{n + 1), and p is an eigenvalue 
of a finite-dimensional matrix (see below). 

The solutions of interest belong to four classes of polynomials, which are often labeled 
/C(s), C{s), A^(s), and Af{s). For every nonnegative polynomial degree n, there are a total 
of 2?T, + 1 number of solutions from these classes. Specifically, defining r = n/2 for n even and 
r = (n — l)/2 for n odd, we have r + 1 solutions in the class /C„(s), n — r solutions in £„(s), 
n — r solutions in A^„(s), and r in A/'„(s). By analogy with spherical harmonics, these 2n + 1 
solutions of degree n are labeled -E^(s), with p (since it is individual to each solution) serving 
double duty as a dummy index from to 2n + 1. Convention holds that these solutions are 
assembled class- by-class in alphabetical order, i.e. for < p < r, E^{s) = JC^is), then the 
next n — r are the solutions C^^is), and so forth. 

For brevity, we describe how to determine only one of the classes of solutions; Romain 
and Jean-Pierre include a complete and detailed description of all four-, and our software 
implementation discussed in Section IIIII can be consulted for full details. A stable algorithm 
introduced by Dobner and Ritter— begins by writing the solution class as 

= €'ns)P^'ns) (21) 

with 

^|J^'^{s) = s--'^ (22) 

and 

Pn^'i^) = Y.h[^-Y^' ■ (23) 
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The pth solution /C^(s) is determined by the coefficients bj in Eq. [231 which represent the 
pth eigenvector of a tridiagonal matrix of the form 



do 90 ■ ■ ■ 
fi di gi ■■■ 
/2 d2 92 ■■■ 



(24) 











dr 



where the nonzero elements fi, di, and gi depend on the solution class, the ellipsoid semi- 
axes, n, and r. 

The above solutions are known as Lame functions of the first kind; they diverge as s — > cxd, 
making it impossible to write solutions to the Laplace equation in exterior regions using these 
functions. A closely related set of solutions do behave appropriately in this limit, and are 
known as Lame functions of the second kind; these only involve the radial-like coordinate A 
and are written 



D. Ellipsoidal Harmonics 

For a given degree n and order p, the interior solid ellipsoidal harmonic is defined as 



(25) 



where 




(26) 




(27) 



the exterior solid ellipsoidal harmonic as 



F^(A,/i,z/) = (2n + l)E^(A,/x,i/)J^(A), 



(28) 



and the surface ellipsoidal harmonic as 




(29) 



The surface harmonics satisfy the orthogonality condition 




'pp' 



(30) 
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where 7^ is a normalization constant; writing the surface integral more explicitly. 

ph rk (,,2 ,,2 

'0 Jh 



ll= I I ' diidu. (31) 



Then the Coulomb potential at r due to a unit charge at r' (with ||r|| > ||r'||) is expanded 
in ellipsoidal harmonics as 



1 -2^1 1 



r — r 



n=0 p=l 



or a similar form if ||r|| < ||r'|| (see, e.g.,-). The normal derivative at the ellipsoid surface 
defined by A = a is computed as^i 

^ A V ^ (33) 

E. Series Solution for Biomolecule Electrostatics 

In analogy with the Kirkwood series solution for a spherical particle^^, we formulate our 
model problem in ellipsoidal coordinates r = (A,/x, v). The potential in region 1, $i(r), can 
be written 

Q 

qk 

k=l 



$i(r) = V +^{r), (34) 

Ci r — Ffc 



where il^{r) is the reaction potential. We expand ip in ellipsoidal harmonics of the first kind 
(interior harmonics) since these are valid within region 1: 

00 2n+l 

^W = EE^nEi:(r). (35) 

n=0 p=l 

Similarly, the potential in region 2 may be expanded in ellipsoidal harmonics of the second 
kind (exterior harmonics), which are regular as ||r|| — > 00: 

00 2n+l 
n=0 p=l 

To determine the constants appearing in these expansions, we apply the boundary conditions 
Eqs. (jl]) and ([5]) by relating them to the Coulomb portion of $1, using the fact that all charges 
are contained inside the ellipsoid (A^ < A). Using the Coulomb-potential expansion from 



Eq.l 

Q Q oo 2n+l 

fc=l ^ I A;=l ^ n=0 p=l 

oo 2n+l 

= E E — ^'(r) (38) 

n=0 p=l 



where 



G'„ = E*^;sEr.(rO. (39) 

fc=l 

Now the first boundary condition, Eq. (jl]), gives us, after equating each (n,p) term in order 
for the relation to hold for all angles, 

^ + B'rMr\ = CI (40) 
GP EP(o) 

— + Bn-EM = C'n (41) 
Ci -TnyCl) 

In order to apply Eq. (I5l), we differentiate each series term by term and equate them, 

'' + -^B^n = C::. (42) 



£2 £2 



dX 1°- 

We can eliminate the coefficients to give the reaction field coefficients in terms of the 
known source charge coefficients G^, 

where, following^, we introduce the notation -E^(A) for logarithmic derivative which respect 
to the argument 

III. IMPLEMENTATION DETAILS 

We have implemented ellipsoidal harmonics in both MATLAB^ and Python. These im- 
plementations are released as open-source software under the Simplified BSD License and 
available online. Computation involves four distinct tasks, which we address in turn be- 
low. We would like to echo Xue and DengS in praising Romain and Jean-Pierre^ for their 
comprehensive discussions and detailed derivations. 
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A. Coordinate Transformations 



Our open source implementations follow a layered design, where high levels use the results 
of lower layers. We begin with the bijection between Cartesian coordinates, and ellipsoidal 
coordinates, from Eq. [T^l We test our implementation by taking a brick of coordinates in 
each octant, convert these coordinates to ellipsoidal coordinates, and then convert back and 
compare to the original Cartesian coordinates. By testing each octant, we assure that we 
have correctly handled the sign ambiguity in Eq. [12] or Eq. [191 

The sign ambiguities associated with squaring the coordinates in these equations can be 
resolved by appealing to the Lame functions of first order, (A, fi, v). These harmonics are 
just the dipoles (e.g. oc a;) , and as such have the same sign as the traditional dipoles. 
Denoting the sign of a coordinate, say x, with the definition Sx = sgn(x), 

Sx = sgn (E?(A)E?(/x)£;°(z/)) (45) 
Sy = sgn{EliX)Elifi)Eliu)) (46) 
s, = sgn{Eli\)E'ME!{i.)), (47) 

however this just moves the sign ambiguity to the Lame functions. We have, from Romain 
Table III or by taking square roots of Eq. [191 



khx = XiiP (48) 
hihy = ^X^ -hl^y? -hl^hl-u^ (49) 



h^kz = y^A2 - hl^hl - f^^^hl - u\ (50) 

where the sign of the square roots must be chosen consistently. We choose the sign of each 
square root as the product of the sign of argument and a sign determined by the semifocal 
axis. Thus we have 

Sx = SxSf,S^ (51) 

Sy = {s\Sh){Sf,Sh){s,ySh) = sxs^SjySh (52) 
Sz = {s\Sk){Sf,Sk){si,Sk) = sxs^s^Sk- (53) 

Multiplying the first and second, and first and third equations, 

SxSy = Sh (54) 

SxSz — Sk- (^^) 
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Now our system has the solution 

sx = SxSySz (56) 
Su = SxSz, (58) 

which we can check using Eq. [5T| 

Sx (^SxSyS;z^(^SxSy^(^SxSz^ (^'^) 

= Sx{SxSy) (60) 

Sz = Sx{SxSz). (61) 

We use the sign assignment 

Sf, = Sh (62) 

Su = Sk- (63) 

when computing the sign of used in the Lame functions, and found in Romain Table II. 
This also gives the simple formula for signs of Cartesian coordinates, 

Sx = sxs^^Su (64) 
Sy = sxSu (65) 
Sz = sxSf,. (66) 

In Python, the EllipsoidalSystem class defines an ellipsoidal coordinate system using 
(a, b, c), and accomplishes the conversion using ellipsoidalCoords () and cartesianCoords () . 
The corresponding functions in MATLAB are approxCartToEll () and ellToCartO. The 
sign of Cartesian coordinates are calculated using Eq. EU and the sign of elhpsoidal coordi- 
nates comes from Eq. [56l There is also a sign ambiguity in the Lame functions, which use 



the sign of \/y? — h? and \fy? — W (our two signs in the derivation above) in the evaluation 
of tl). Thus, the calcLameO, also calcLameO in MATLAB, method also takes the sign of 
[i and V, regardless of the coordinate used for evaluation. 



B. Evaluating the First-Kind Lame Functions 



The evaluation of Lame functions using is straightforward, following Romain; in the 
Python implementation, see the calcLameO method. The (n,p) is first converted to a 
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solution classs (/C, £, A^, or M), and a class-specific index p' into the solutions in that class. 
For instance, for n = 2 there exist 5 solutions: 2 IC2 and one of each of the others, so (2, 2) 
is associated with £2- 

The factor can then be calculated trivially; Romain Table 11 gives the factors for all 
four solution classes^. Next, we assemble the tridiagonal matrix for the given Lame type 
and solve the eigenproblem. Since the matrix is well scaled, the problem is easily solved 
using the default LAPACK routine. We use the p'th eigenvector to provide the coefficients 
for the polynomial of Eq. [221 Following standard practice, we normalize the eigenvector so 
that the highest power of the argument A in the sum has coefficient unity; this is prescribed 
in Dassios^^ and elsewhere. Derivatives of Lame functions simply invoke the product rule 
for Eq. [211 and derivatives of the appropriate ip and polynomial P are simple to compute. 

C. Evaluating Integrals for Second-Kind Lame Functions 

Lame functions of the second kind, given by Eq. [251 require the calculation of from 
Eq. |26l Although this quantity can in principle be calculated from elementary elliptic 
integrals^, the decomposition is ill-conditioned and thus we choose to evaluate the integral 
numerically. Adaptive quadrature rules, such as quadO from Matlab, are quite efficient for 
this job. Since the argument appears as the upper limit of integration, the derivative of 
amounts to a single evaluation of the integrand. 

D. Computing the Normalization Constants 

Apparently, the greatest numerical challenge is presented by evaluation of the normal- 
ization constant 7^ used for surface ellipsoidal harmonics and representing the Coulomb 
potential. Romain and Jean-Pierre present a way to evaluate these constants using four 
basic elliptic integrals. Each integral has an integrable singularity at the end point, and 
thus requires something like a Gauss-Kronrod quadrature in order to evaluate the integral 
accurately through the use of transformations. The simple midpoint quadrature currently 
used in the Python implementation converges quite slowly. Both implementations have been 
checked against analytic answers provided by Dassios^i for n = 0, 1, 2; these low-order terms 
are correct to at least five significant figures. 
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A further complication arises in computing the expansion for the fundamental solution 
of Laplace's equation in Eq. [321 For certain geometries (depending on the positions of the 
points and the ellipsoid semi-axes), the value of can become very large, and is balanced 
by a very small 7^. Both quantities are computed numerically, and so this cancellation 
can quickly render the computation inaccurate. As a result, expansions beyond order 11-12 
are often computed to unsatisfactory accuracy (see Fig. [1]), although the precise breakdown 
depends on the semiaxes and position of the evaluation point. It seems clear that this 
computation should be reformulated in order to calculate the quotient directly. 

IV. RESULTS 

We have verified the two implementations against one another, and validated the code 
using the extensive list of identities presented by Dassios and Kariotoi>2i, which were in- 
valuable during development and testing. In this section, we present three sets of results to 
illustrate that the calculations are indeed correct. For simplicity in exposition, and due to 
the Python numerical integration challenges discussed above, results are from the MATLAB 
implement at ion . 

First, we illustrate that the series expansion of the Coulomb potential converges to the 
exact result with increasing order, i.e. that the error associated with truncating the infinite 
sum in Eq. |32] goes to zero. Figure [1] plots the magnitude of this error in the case of 
ellipsoid with semi-axes a = 2, b = 1.5, c = 1, with internal charge r = (0,0,0.5) and field 
point r' = (0, 0, 2); three of the sub-figures illustrate the challenge of accurate computation, 
plotting the magnitudes of E^, F^, and 7^. Numerical error in the required integrations 
entails a limited accuracy — in this case, a relative error of approximately 10""^. Other test 
examples exhibit slower convergence and more complex non-monotonicity. 

Our second result demonstrates that the electrostatic solvation free energy of a single 
charge at the origin recovers the well-known result for the sphere limit, also known as the 
Born energy, as we model an ellipsoid with semi-axes a = 1 + A, b = l + A/5, c = 1 + A/lO 
(all units in Angstrom), and let A approach zero (Figure [2]). These calculations show that 
the implementations behave correctly for small but finite differences between the ellipsoid 
semi-axes. Unfortunately, vanishingly small differences pose numerical challenges (note that 
the ranges for /x and u depend on the differences in the semi-axes). Obtaining the spherical 
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harmonics in the actual hmit requires careful analysisi*^. 

Finally, using the same dielectric constants as in the previous example, we calculate the 
electrostatic solvation free energy of a single charge located at (3, 4, 5) inside a protein-like 
ellipsoid with semi-axes a = 15, b = 12, and c = 10. For this system, we can validate 
the implementation numerically using simple BEM calculations, which is also written in 
MATLAB and included with the ellipsoidal harmonics software. Figure [3] plots the deviation 
between the semi-analytical result and the numerical BEM results, as a function of the 
number of boundary elements (flat triangles) used to approximate the ellipsoid surface. 
We observe the expected linear convergence to the semi-analytical result, which suggests 
the correctness of our ellipsoidal-harmonic calculation. In the present implementation, more 
accurate validation is impractical because meshes are generated from within MATLAB using 
the built-in function ellipsoid, and BEM calculations are performed using dense 0{N'^) 
algorithms. Fast BEM solvers^i^ will be needed for more stringent testing. Also, as Fig. [1] 
illustrates, more accurate methods need to be found to compute the needed integrals. 



V. DISCUSSION 



In this paper we have presented two open-source implementations of ellipsoidal harmon- 
ics to allow rapid, analytical calculations of many problems of potential theory, including 
electrostatics, electromagnetics, elasticity, and fluid mechanics. The implementations are 
written in MATLAB and Python, and released under the Simplified BSD License, and are 
freely available online ??. These development efforts stemmed from our interest in im- 
proving continuum Poisson-based models of the electrostatic interactions between biological 
molecules and the surrounding solvent^i^i^. Many important studies have relied on analyt- 
ically solvable spherical geometries^, and a variety of state-of-the-art approximate models 
rely on analyses in spherical harmonics, e.g.— i^. We hope that the present implementation 
of ellipsoidal harmonics, or at least our description of pitfalls we encountered along the way, 
may provide a starting point for more accurate Generalized-Born models^"^ or other fast 
approximations to the Poisson problent^""— . 

We do not claim to have covered the entire field of ellipsoidal harmonics, however, or 
that our implementations are production-level codes with all mathematical subtleties fully 
addressed. The implementations released with this paper remain under active development, 
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and interested readers are encouraged to contact us with questions, requests, bug reports, 
and suggestions of all kinds. Currently, the most important area for improvement appears 
to be in handling coordinate transformations, particularly methods for accurately mapping 
Cartesian to ellipsoidal coordinates'', and transformations that seem to need negative A where 
other publications indicate that non-negative A is adequate. Our results also illustrate the 
need for more accurate ways to evaluate integrals for the second-kind Lame functions and 
the normalization constants. 

The present implementation includes the basic algorithms needed to compute the ellip- 
soidal harmonics, and we have used these primitives to address a fairly simple problem of 
potential theory (a mixed-dielectric Poisson problem). The primitives can be used for a 
variety of other problems in potential theory, which are of interest on scales ranging from 
atoms^ii^ to entire planets and the solar system'^'^^. In addition, recent work by Ritter 
and colleagues has elucidated the relationship between ellipsoidal harmonics and boundary- 
integral operators of potential theory for ellipsoidal boundaries (e.g.,— i^). Our approach 
to approximating biomolecule electrostatic problems using boundary-integral methods^"— 
should benefit significantly. The most recent analysis employed spherical harmonics and 
found that accurate approximations could be obtained using parameters that correspond 
not to spheres but surfaces close to spheres^. We conducted this work because ellipsoids 
have been used for a number of other molecular electrostatic models^"— and thus the 
boundary-integral analysis of Ritter et al. may allow new approach to approximating the 
Poisson problem. Also, much of the earlier work employed only the lowest-order modes 
(which can be expressed in closed form), and did not allow extension to the higher-order 
ones that must be computed. 

In closing: ellipsoidal harmonics have a wide range of applicability, but historically the 
mechanical details required to use them have been an unfortunate deterrent. Modern com- 
puters and numerical algorithms can easily handle these details to enable powerful new 
modeling tools, and in this work we take a first step towards this goal. It is already clear 
that computing these functions can motivate more fundamental numerical research, and we 
hope others will find this subject and its applications equally fascinating. It seems entirely 
likely that mechanical computation derailed the development of other powerful formalisms in 
the era before digital computers; perhaps computational science should undertake a journey 
of re-discovery. 
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FIG. 1. Convergence of the ellipsoidal-harmonic approximation of the Coulomb potential as a 
function of the degree n; the ellipsoid is defined hy a = 2, b = 1.5, c = 1, the source charge is at 
(0,0,0.5), and the field point is (0,0,2). 
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FIG. 2. Validating the semi-analytical calculation of electrostatic solvation free energies using the 
analytical result for a sphere with a central charge (a quantity known as the Born energy). We 
define A > and an ellipsoid with semi-axes a = l-|-A, 6=1-|- A/5, and c = 1 + A/10. A 
single +le charge is situated at the origin, and we have ei = 4 and e2 = 80. As A — )■ 0, the 
semi-analytical ellipsoidal results converge to the exact energy for the sphere, which supports the 
correctness of our implementation. 
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FIG. 3. Validating the semi-analytical calculation of electrostatic solvation free energies using 
numerical simulations based on the boundary-element method (BEM). The ellipsoid has semi-axes 
a = 15, 6 = 12, and c = 10, and a single +le charge is situated inside at (3,4,5); as in Fig. [2l 
ei = 4 and €2 = 80. The BEM results converge linearly as a function of the number of unknowns 
(one per panel), which indicates that the semi- analytical method is returning the correct result. 
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